Non-Universality of the Specific Heat in Glass Forming Systems 
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We present new simulation results for the specific heat in a classical model of a binary mixture 
glass-former in two dimensions. We show that in addition to the formerly observed specific heat 
peak there is a second peak at lower temperatures which was not observable in earlier simulations. 
This is a surprise, as most texts on the glass transition expect a single specific heat peak. We explain 
the physics of the two specific heat peaks by the micro-melting of two types of clusters. While this 
physics is easily accessible, the consequences are that one should not expect any universality in the 
temperature dependence of the specific heat in glass formers. 



The thermodynamic properties of glass-formers near 
the glass transition have been a subject of intensive and 
far from settled research for more than half a century [1- 
4] . The temperature dependence of the entropy and how 
the entropy extrapolates to low temperatures gave rise to 
the so-called Kauzmann paradox [1] that remains confus- 
ing to the present time. Important to the understanding 
of these issues is the specific heat, either at a constant 
volume or at a constant pressure, since its integral over 
a temperature path provides the entropy. Experimental 
measurements of the specific heat in glass-forming sys- 
tems are obtained from the linear response to either slow 
cooling (or heating) or to oscillatory perturbations with 
a given frequency about a constant temperature. The 
latter method gives rise to a complex specific heat with 
the constraint that the zero frequency limit of the real 
part should be identified with thermodynamic measure- 
ments. Such measurements reveal anomalies in the tem- 
perature dependence of the specific heat, including the 
so called "specific heat peak" in the vicinity of the glass 
transition. In fact, throughout the literature on the glass 
transition one finds references to the specific heat peak 
[5]. In this Letter we show that this concept must be 
discarded, since depending on the detailed physics of the 
system there can be two or multiple specific heat peaks. 
We will present evidence for a model system with two 
specific heat peaks, explain in detail the physical origin 
of the latter, and point out the important consequence 
that there is very little (or no) universality that can be 
expected in the thermodynamic properties of different 
glass-formers. 

The model discussed here is the classical example [6, 7] 
of a glass-forming binary mixture of N particles in a 2- 
dimensional domain of area V, interacting via a soft 1 / r 12 
repulsion with a 'diameter' ratio of 1.4. We refer the 
reader to the extensive work done on this system [6-10]. 
The sum-up of this work is that the model is a bona fide 
glass-forming liquid meeting all the criteria of a glass 
transition. 

In short, the system consists of an cquimolar mixture 
of two types of point-particles, "large" with interaction 
range a% — 1.4 and "small" with interaction range o~\ = 
1, respectively, but with the same mass m. In general, 
the three pairwise additive interactions are given by the 



purely repulsive soft-core potentials 



<f>ab(r) = e (j^r) 



a, b= 1,2 



(1) 



where a aa = a a and a a b = (o~a + c&)/2. The cutoff radii 
of the interaction are set at 4.5cr a (,. The units of mass, 
length, time and temperature are m, <T\, t = ai^/m/e 
and e/fcs, respectively, with ks being Boltzmann's con- 
stant. In numerical calculations the stiffness parameter 
of the potential (1) was chosen to be n = 12. 

The isochoric specific heat is determined by the fluctu- 
ations of the energy of the system at a given temperature: 
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The specific heat of our binary mixture model was mea- 
sured at constant volume in [7, 11] and by us. We have 
used the last equation which allows one to estimate the 
specific heat in a single run of the canonical ensemble 
Monte Carlo simulations. At each temperature the den- 
sity was chosen in accordance with the simulation results 
in an NPT ensemble as described in [7] with the pres- 
sure value fixed at P — 13.5. As the initial configuration 
in the Monte Carlo process the last configuration of the 
molecular dynamics run for this model at given temper- 
ature after 1.3 x 10 8 time steps was used. After short 
equilibration the potential energy distribution functions 
were measured during 2 x 10 6 Monte Carlo sweeps. The 
acceptance rate was chosen to be 30%. Simulations were 
performed with N = 1024 particles in a square cell with 
periodic boundery conditions. Our simulations appear to 
provide trustable values of Cy down to lowest tempera- 
tures where the value of the specific heat coincides with 
that of two-dimensional solid, i.e. Cy = 2. What could 
not be seen in earlier simulations is that there is a smaller 
second peak of the specific heat at lower temperatures. 
To resolve it to the naked eye we present in Fig. 1 a 
blow-up of the region of lowest temperatures where the 
second peak is more obvious. 

To understand the nature of the specific heat anomalies 
one must understand the physics that is behind the glassy 
behavior of this model in general and the existence of the 
two specific heat peaks in particular. When the tempera- 
ture is lowered at a fixed pressure this system [12] (as well 



FIG. 1: Color online: in dots: the temperature depenc 
of the specific heat in the binary mixture model, comr. 
at constant volume, such that the volume agrees with the 
pressure P=13.5 [7] at each temperature. The data indicate 
the existence of two specific heat peaks, one prominent at 
about T = 0.5 and a smaller on at about T = 0.1, and see the 
inset for finer detail. The data at lower temperature represent 
two years of computing time and are believed to be trustable. 



as many other glass-formers [13-18] tends to form micro- 
clusters of local order. In the present case large particles 
form hexagonal ordering first (starting at about T = 0.5, 
and at lower temperatures (around T = 0.1) also the 
small particles form hexagonal clusters. The clusters are 
not that huge, with at most 0(100) particles, (cf. Fig. 2), 
depending on the temperature and the aging time. But 
we have shown that the long time properties of correla- 
tion functions are entirely carried by the micro-clusters 
[12]. Below we will refer to the micro-clusters as curds 
and the liquid phase as whey. We will argue that the 
specific heat responds to the micro-melting of the clus- 
ters - those of small particles at the lowest temperatures 
and those of the larger particles at higher temperatures. 
The large increase in the number of degrees of freedom 
when a particle leaves a crystalline cluster and joins the 
liquid background is the basic reason for the increase in 
entropy that is seen as a specific heat peak. 

We will calculate the specific heat at constant volume 
per particle from the exact expression that can be derived 
[19] for any system with inverse power law potential r~ n , 
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For our system n = 12, K is the bulk modulus and K c 
[20, 21] is given by: 
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The bulk modulus requires an equation of state for its 
calculation. In the rest of this Letter we will therefore 



FIG. 2: (Color online). A snapshot of the system at T = 0.44. 
In colours we highlight the clusters of large particles in local 
hexagonal order. The colours have no meaning. 



derive an approximate equation of state and compute the 
specific heat, exposing the origins of the two peaks. 

To start we define v l w , v^,, v l c and vl respectively as 
the volume of large particle in the whey, small particle in 
the whey, large particle in the solid and small particle in 
the solid. Similarly we denote by e l w , e^, e l c and e s c the 
energy of a large and small particle in the in the whey 
and in the crystalline phase respectively. Needless to say, 
all these quantities are temperature and pressure depen- 
dent; we will therefore explicitly use our low temperature 
knowledge concerning v c and v s c in the crystalline phase, 
but treat the difference v l w — v l c and v s w — v s c as constants 
that we estimate below from our simulation knowledge. 
Similarly we estimate e l c and e s c from our knowledge of the 
hexagonal lattices at T = 0. We assume that e l w w and 
similarly ~ e s c since our simulations indicate a very 
small change in these parameters, see Table I. It should 
be stressed that the enthalpy change at these pressures 
are almost all due to the PV term. This will result in a 
semi-quantitative theory ascribing the important changes 
in specific heat to the changes in the fraction of particles 
in curds and whey. In other words the number of parti- 
cles in the whey and the number of clusters are all explicit 
functions of temperature and pressure. 
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3.69 


2.07 


3.76 


2.16 


1.43 


0.92 


1.58 


0.94 



TABLE I: Parameters used in the calculation of the specific 
heat 

As the condensed phase consists of clusters of large and 
small particles, we use the notation for the number of 
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clusters of n large particles and for the clusters of m 
small particles. Here we only need the intensive variables 
Pi = 2J2 n Ni/N, pi = 2j2 m Ni/N pi = 2Ni/N and 
p b w = 2N^/N which stand for the fraction of large parti- 
cles and small particles in the curds, and large particles 
and small particles in the whey, such that pr c + p^ — 1 
and p s c + p s w = 1 . Using these variables we can write an 
expression for the volume per particle v = V/N: 



-Pi 



-Pc 



(5) 



At this point we need to derive expressions for p l c and 
p s c . To do so we need to remember that in the rel- 
evant range of temperatures the large particles in the 
whey can occupy either hexagonal or hcptagonal Voronoi 
cells, whereas small particles can occupy only pentago- 
nal or hexagonal cells [8, 9, 12]. Accordingly there are 
<7~ « (2 6 — l)/6 + 2 7 /7 ways to organize the neighbours 
of a large particle in the whey (neglecting the rare large 
particle in heptagonal neighbourhood), but only one way 
in the cluster. Similarly, there are g^ « (2 6 - 1)/6+ 2 5 /5 
ways to organize a small particle in the whey. We note 
that this estimate assumes that the relative occurrence 
of the different Voronoi cells is temperature independent. 
While reasonable at higher temperatures [12], at lower 
temperature one should use the full statistical mechanics 
as presented in [9] to get more accurate estimates of g l w 
and g!^. This is not our aim here; we aim at a physical 
understanding of the specific heat peaks rather than an 
accurate theory. We thus end up with the simple esti- 
mates 
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It is important to note that the combination of Eq. 
(5) together with Eqs. (6) and (7) provides a mechanical 
equation of state. We will now compute C v directly from 
Eq. (3). The peaks in the specific heat are determined 
by the temperature dependence of p e c (P, T) and p s c {P 1 T) 
each which has a temperature and pressure derivatives 
that peaks at a different temperature, denoted as T e (P) 
and T S (P). As said above we take Av l = v l w — v c and 
Av s = — v s c as approximately constants (as a function 
of temperature and pressure). The constants are esti- 
mated from the condition that the second temperature 
derivative of p^,(P,T) and p s c {P,T) should vanish. From 
this conditions we find 



At 



T*(P*)lngi/P* , Av s « T*(P*) In <fJP* 



(8) 

where P* is the pressure for which the peaks in the 
derivatives are observed (13.5 in our simulations). This 
is equivalent to a linear dependence of the specific heat 
peaks as a function of pressure, T l {P)/T\P*) = P/P* 
and similarly for the small particles. 
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FIG. 3: Specific heat at constant volume as predicted by the 
simple theory which is based on the mechanical equation of 
state supplied by Eqs. (5) and (6) and (7). Note that the 
theory predicts the two peaks which are associated with the 
micro-melting or micro-freezing of the clusters of large and 
small particles respectively. The magnitude of the peaks is 
too high, reflecting terms missing in the simple approach, 
like the effect of anharmonicity at the lowest temperatures 
which are negative, tending to decrease the height of the low- 
temperature peak. 



In terms of these objects we can rewrite 



v = v c (P,T)+Av t (l 
dv\ ( dv r 
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To compute the temperature dependence of (§p) T we 
need first to determine its T — > limit, which is deter- 
mined by the first term on the RHS of Eq. (9) as the 
other terms on the RHS decay exponentially fast when 
T — > 0. Since we have already exact results for the bulk 
modulus for the present model, we return to Eqs. (3) 
and (4). We know on the one hand that limT^o C v = 2 
and that (U)/N w 2.94 over the whole interesting tem- 
prature range, cf. [19]. The compressibility n is related 
to the bulk modulus via k — — (-§p) T /v = 1/K and 

therefore {^p) T ~ —1/(123.5 — 35T) is easily estimated 
as T — y 0. Eor simplicity we will use this approximation 
up to T w 0.5. 

Having all the ingredients we can compute C v /N . The 
parameters used were estimated from the numerical sim- 
ulation and are summarized in Table I. Since the aim of 
this subsection is only semi-quantitative, we do not make 
any attempt of parameter fitting, and show the result of 
the calculation in Fig. 3. 

Indeed, the theoretical calculation exhibits the exis- 
tence of two, rather than one, specific heat peaks. We 
can now explain the origin of the peaks as resulting from 

the derivatives (§f) an d ■ These derivatives 

change most abruptly when the micro-clusters form (or 
dissolve), each at a specific temperature determined by 
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{h s w —h s c )/\ng s w and (h^ — h^/lng^. Note that there can 
be pressures (both upper and lower boundaries) where 
the the sign of {h s w — h s c ) or {h l w — h e c ) change sign and 
the peak can be lost. 

In summary, we have presented a simulational discov- 
ery of an unexpected second peak in the temperature 
dependence of the specific heat of a popular model of 
glass-formation. This discovery means that the universal 
expectation of seeing a single "specific heat peak" should 
be seriously revised. As there are two peaks in this ex- 
amples, other examples may have multiple peaks or one 
peak. Even in this system at a different pressure the 
two peaks may merge, giving the appearance of a single 
peak. We have also presented a detailed explanation of 
the existence of the two peaks, based on an approximate 



equation of state that we have derived on the basis of sta- 
tistical mechanics. Indeed, it is not difficult to explain 
the two peaks in terms of the effects of micro-melting 
of clusters of small and than of large particles when the 
temperature is increased. The important consequence is, 
however, that no universality is expected for the thermo- 
dynamics properties of different glass formers. 
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